clear all
load env23_12.mat

SNR = 1:0.25:6;
contBGW = zeros(length(SNR),k+2);

for snr = 1:length(SNR)
    
    cont = 0;
    rand('state',0);
	randn('state',0);
    
    tx = [];
    rxCI = [];
    rxMLD = [];
    
    while (cont < 1e4)
    
        % Mensagem
        u = randint(1, k);
    
        % Palavra Binaria Codificada
        cw = abs(mod((u * G), 2));
        
        % Acumula mensgens transmitidas
        tx = [tx cw];

        % Modulação BPSK
        x = 2*cw - 1;

        % Ruido AWGN
        y = awgn(x, SNR(snr));
        
        % Zera distâncias Euclidianas
        dEuc = [];
        
        %------------------------------------------------------------------
        % CI
        %------------------------------------------------------------------
        % Inicia com primeiro elemento de select_list
        line = 1;
        
        % Ordena itens dos mais confiáveis para os menos confiáveis
        [val, ordem] = sort(abs(y),'descend');
        
        % Determinante de Gparcial deve ser ímpar
        detGp = 0;
        while (mod(detGp, 2) == 0)
        
            for i = 1:k
                k_bits(i) = ordem(1, select_list(line,i));
            end
            
    
            for i = 1:k
                msg0(1, i) = y(1,k_bits(1,i))>0;
                Gparcial(:,i) = G(:,k_bits(1,i));
            end
            
            detGp = det(Gparcial);
            line = line + 1;
            
            if (line > 41)
                disp('Impossível gerar matriz Gnova!');
            end
        end
        
        % Calcula Gnova
        Gnova = abs(mod((round(det(Gparcial)*inv(Gparcial)*G)),2));
        [r, c] = size(Gnova);

        masc = [zeros(1,k-1) 1];
        stopCond = 0;
        
        % As outras candidatas são permutações de cada bit da palavra
        % original
        for j = 1:(k+1)
            
            if (j == 1)
                % Primeira palavra candidata
                msg_temp = msg0;
            else
                % Demais palavras candidatas
                msg_temp = xor(msg0, masc);
                masc = circshift(masc',-1)';
            end
            
            cw_temp = abs(mod((msg_temp * Gnova),2));

            cw_candidatas(j,:) = cw_temp;
            dEuc(j) = sum((y - (2*cw_temp-1)).^2);
                    
            % Limitante de BGW
            somador = (-2*cw_temp + 1).*y;
            somador = sort(somador,'descend');        
            limBGW = sum(somador(1:dHmin));
            if (limBGW < 0)
                stopCond = 1;
                decodCI = cw_temp;
                contBGW(snr,j) = contBGW(snr,j) + 1;
                break;
            end

        end
        
        % Se limitande BGW não atua, escolhe a palavra com minima distancia
        % Euclidiana
        if (stopCond == 0)
            % Escolhe-se a palavra com menor distancia encludiana de y
            [val, pos] = min(dEuc);
            decodCI = cw_candidatas(pos,:);
            contBGW(snr,14) = contBGW(snr,14) + 1;
        end
    
        % Acumula mensagens recebidas
        rxCI = [rxCI decodCI];
        
        
        %------------------------------------------------------------------
        % MLD
        %------------------------------------------------------------------
        % Procura MLD
        [r,c] = size(palavrasCod);

        % Calcula-se uma MLD inicial, para o primeiro elemento de palavrasCod
        pc_elem = palavrasCod(1,:);
        %mld = sum(xor(pc_elem, d));
        mld = sum((y - (2*pc_elem-1)).^2);

        % A palavra decodificada assume esse elemento inicial
        decodMLD = pc_elem;

        for i = 2:r
            % Para cada linha de palavrasCod, calcula-se a distancia
            % Euclidiana
            pc_elem = palavrasCod(i,:);
            dEucMLD = sum((y - (2*pc_elem-1)).^2);
    
            % Se a dHam calculada for menor do que a MLD calculada, seleciona-se a
            % palavra codigo correspondente como palavra decodificada
            if (dEucMLD < mld)
                decodMLD = pc_elem;
                mld = dEucMLD;
            end
        end
        
        
        % Acumula mensagens recebidas
        rxMLD = [rxMLD decodMLD>0];

        cont = cont + 1;
    end
    
    ber23_12ci(snr) = sum(xor(tx, rxCI))/(cont*n)
    ber23_12mld(snr) = sum(xor(tx, rxMLD))/(cont*n)
end

save simul23_12.mat ber23_12ci ber23_12mld SNR contBGW

% Resultados
figure(1);
semilogy(SNR, ber23_12mld,'b');
hold all;
semilogy(SNR, ber23_12ci,'r');
hold off;
title('Golay code (23,12)');
xlabel('SNR (dB)');
ylabel('BER');
legend('MLD','CI');
axis([min(SNR) max(SNR) 1e-6 1e-1])
grid on;